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We develop an efficient method to calculate probabilities of large deviations from the typical be- 
havior (rare events) in reaction-diffusion systems. The method is based on a semiclassical treatment 
of underlying "quantum" Hamiltonian, encoding the system's evolution. To this end we formulate 
corresponding canonical dynamical system and investigate its phase portrait. The method is pre- 
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I. INTRODUCTION 

Reaction-diffusion models have a vast area of 
applications 1-2 ranging from kinetics of chemical 
reactions 3 , biological populations 4 **^ and epidemics^ to 
the dynamics of financial markets** and ecology^. The 
models describe dynamics of a number of particles whose 
reactions are specified by a certain set or rules. The rules 
have a probabilistic nature and are most conveniently for- 
mulated on a lattice in a d-dimensional space. We shall 
restrict our attention to a wide subclass of such models, 
where the particles execute random walks (diffuse) on the 
lattice, while the reactions between them are purely lo- 
cal (on-site). Once the lattice, reaction rules and initial 
conditions are specified, one is interested to find statisti- 
cal characteristics of the system's subsequent evolution. 
This goal may be accomplished with various degrees of 
detailing and accuracy. 

The most detailed information is contained in the prob- 
ability distribution functions (PDF) of every possible mi- 
croscopic state of the system. The PDF is a solution of 
an exponentially large system of Master equations, which 
specify probabilities of transition between every two mi- 
croscopic states of the system. Analytical solution of 
the Master equations is usually unrealistic and besides 
the information contained in them is excessive. There- 
fore various approximation schemes are in order. The 
simplest one is the mean-field approximation, where a 
closed set of equations for average quantities (e.g. con- 
centrations) is obtained by an approximate decoupling of 
higher moments. The mean-field theory describes a typ- 
ical evolution of the system, if the fluctuations are weak 
in a certain sense. Probability of small deviations from 
the mean-field predictions may be found with the help 
of Fokker-Planck (FP) equation, ft substitutes the dis- 
creet Master equation by a continuum (biased) diffusion 
equation in the space of concentrations. Analysis of FP 
equation is usually complicated*^, moreover the approx- 
imation is reliable for small deviations only and fails to 
provide probability of large deviations from the typical 
evolution. 

Much attention was attracted recently to reaction- 
diffusion systems that are in a close proximity to dy- 
namic phase-transitionsii*i£ii£ (for recent reviews see 



e.g. Refs. [TUl). By fine-tuning one of the parameters 
some systems may be brought to a point of quantitative 
change of their behavior (e.g. stable finite concentration 
vs. extinction). In a vicinity of the transition, neither 
mean-field nor FP can accurately predict the long-time 
scaling of the system's characteristics, such as e.g. parti- 
cles concentration. The field-theoretical renormalization 
group (RG) methods were developed and successfully ap- 
plied to a number of examplesiSiiLiSiiS. In particular, 
the directed percolation universality class was identified 
and studied as the most robust universality class for the 
dynamic phase transitions 14-20 ' 21,22 . 

In the present work we address somewhat different set 
of questions. We consider a generic reaction-diffusion 
system that either does not exhibit, or is far enough from 
the phase transition. A typical evolution scenario and 
probability of small deviations are well described by the 
mean-field theory and the FP equation. We shall look, 
however, for a probability of large deviations from the 
typical behavior. "Large" deviation may be loosely char- 
acterized as being of the same order (or larger) as the 
typical value (as opposed to the "small" one, which is of 
the order of the square root of the typical value) . Since 
the occurrence of such large deviations has a very small 
probability, they my be dubbed as "rare events" . Despite 
being rare the "rare events" may be of great interest, es- 
pecially if they cause extreme consequences. Some of the 
examples include: proliferation of virus after immuniza- 
tion (causing death of a patient); large fluctuations of 
number of neutrons in a nuclear reactor (causing explo- 
sion), etc. Clearly in these and many other examples one 
is interested to know rather precisely how improbable are 
improbable events. 

Here we develop a rigorous, simple and efficient method 
to calculate the rare events statistics in reaction-diffusion 
systems. To this end we develop a Hamiltonian formula- 
tion of reaction-diffusion dynamics. Although the system 
is specified by a set of rules, rather than a Hamiltonian, 
one may nevertheless show that there is a certain canon- 
ical Hamiltonian associated with the system's dynamics. 
More precisely, the Master equation may be reformu- 
lated as "quantum" (many-body) Schrodinger equation 
with some "quantum" Hamiltonian. This observation is 
not new and is sometimes referred to as Doi's operator 
techniqu e - 23 ! 24 . In fact, its "quantum" version is the basis 
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for the field-theoretical RG treatment of the dynamical 
phase transitionsi£iiLi£ii&. Here we notice that the clas- 
sical (or rather semiclassical) dynamics of the very same 
Hamiltonian carries a lot of useful information about 
reaction-diffusion systems. In particular, it provides all 
the information about the rare events statistics. To ex- 
tract this information, it is convenient to formulate the 
underlying Hamiltonian in classical terms (as function of 
momenta and coordinates) , rather than creation and an- 
nihilation operators, as is customary in the "quantum" 
approach 23 i 24 i 25 . 

A particularly convenient tool to visualize the sys- 
tem's dynamics is a phase portrait of the corresponding 
Hamiltonian. It consists of lines (or surfaces) of constant 
"energy" (the integral of motion naturally existing in a 
Hamiltonian system) in the space of canonical momenta 
and coordinates. The mean-field (typical) evolution cor- 
responds to a particular manifold of zero energy, given by 
fixed value of the canonical momenta, p = 1. Rare events 
may be specified by certain initial and finite conditions in 
the phase space of the dynamical system, which, in gen- 
eral, do not belong to the mean-field manifold. Probabil- 
ity of the rare event is proportional to exp{— S 1 }, where 
S is the classical action on a unique trajectory, satisfy- 
ing the specified boundary conditions. The problem is 
therefore reduced to finding an evolution of the classical 
dynamical system, whose quantized Hamiltonian encodes 
the Master equation. Such task is substantially simpler 
than solving the full "quantum" Master equation. In 
fact, even probability of small deviations is much more 
efficiently calculated in our semiclassical method, than 
via solution of the FP equation (though the latter is also 
applicable). For large deviations, however, the FP ap- 
proach leads to inaccurate results, while the semiclassi- 
cal method provides the simple and accurate prescrip- 
tion. The very similar strategy was recently applied for 
the calculation of the full current statistics of mesoscopic 

In this paper we develop the semiclassical method us- 
ing a number of reaction-diffusion models as examples. 
We tried to keep the presentation self-contained and ped- 
agogical. We start in section [H] from the model of binary 
annihilation in zero dimensions. In section lTTTl wc compli- 
cate the model by including branching and discuss extinc- 
tion probability of a system having a stable population in 
the mean-field approximation. Section llVI is devoted to 
the extension of the formalism to a d-dimensional space. 
As an example we find an extinction probability of a fi- 
nite cluster. In section a population dynamics model 
with three reaction channels: reproduction, death and 
emigration is considered in a d-dimensional space. The 
model possesses a long lasting meta-stable state with a 
fixed population, that eventually escapes into the state 
of unlimited population growth. We show how the semi- 
classical method may be used to calculate the lifetime 
of such meta-stable state. Finally some conclusions and 
some open problems are discussed in section IVll 



II. BINARY ANNIHILATION 

The simplest reaction, which we use to introduce no- 
tations and set the stage for farther discussions, is the 
binary annihilation process. It describes a chemical re- 
action, where two identical particles, A, form a stable 
aggregate with the probability A whuch does not involve 

in further reactions: A + A — > 0. We start from the 
zero-dimensional version of the model, where every par- 
ticle may react with every other. Such reaction is fully 
described by the following Master equation: 

j t Pn{t) = ~ [(n+ 2)(n + l)P n+2 (t) - n(n - l)P n (t)} , 

where P n (t) is a probability to find n particles at time 
t. The Master equation is to be supplemented with an 
initial distribution, e.g. P„(0) = e"™°njf/ji! - the Poisson 
distribution with the mean value no, or P„(0) = S n ^ l0 - 
the fixed initial particle number. Let us define now the 
generating function as 



(2) 



Knowing the generating function, one may find a proba- 
bility of having (integer) n particles at time t as P n (t) = 
dpG(p 7 t)\ p=0 /nl. If n 3> 1 it is more convenient, to use 
an alternative representation: 



(3) 



where the integration is performed over a closed contour 
on the complex p - plane, encircling p = and going 
through the region of analyticity of G(p,t). 

The point p = 1 plays a special role in this formulation. 
First of all, the conservation of probability demands the 
fundamental normalization condition: 



G(l,t) = 1 



(4) 



Second, the moments of the PDF, P n (t), may be ex- 
pressed through derivatives of the generating function at 
p=l, e.g. (n(t)) =EnPn(t) = d p G(p,t)\ p=1 . 

n 

In terms of the generating function the Master equa- 
tion may be identically rewritten as 



dG A. o <9 2 G 



(5) 



This equation is to be solved with some initial condition, 
e.g. G(p, 0) = exp{rt (p — 1)} for the Poisson initial dis- 
tribution or G(p, 0) = p n ° for rigidly fixed initial particle 
number. The solution should satisfy the normalization 
condition, Eq (0J, at any time. In addition, all physi- 
cally acceptable solutions must have all p-derivatives at 
p = non-negative. 
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One may consider Eq. 10 as the " Schrodinger" equa- 
tion: 



dt 



G = -HG , 



(0) 



where the "quantum" Hamiltonian operator, H, in the p 
("momentum") representation is: 



(7) 



Here we have introduced the "coordinate" operator, q as 



q 



d_ 

dp ' 



\p,q] = l 



(8) 



The "Hamiltonian", Eq. (Q, is normally ordered and not 
Hermitian. The last fact does not present any significant 
difficulties, however. 

If the "quantum" fluctuations are weak (which in 
present case is true as long as (n(t)) ^> 1), one 
may employ the WKB approximation to solve the 
" Schrodinger" -Master equation. Using anzatz G(p,t) = 
exp{— S(jp, t)} and expanding S(p, t) to the leading order 
in 1/A, one obtains the classical Hamilton-Jacoby equa- 
tion: 



dS_ 



H[p, 



OS 
dp 



A 



:(P 2 



1) 



as 

dp 



(9) 



Instead of directly solving the Hamilton-Jacoby equation 
we'll develop the Hamilton approach, which is much more 
convenient for finite dimensional applications. 

To this end we employ the Feynman path integral rep- 
resentation, which may be derived introducing resolution 
of unity at each infinitesimal time-step and employing 
the normal ordering. As a result, one finds for the gen- 
erating function: 



G(p,t) 



lim 

M^oo 



n' dp k dq k 
—■ e 



S[pk,qk] 



k=o 



2tt 



(10) 



where the discrete representation for the action S[pk, qk] 
is given by 



S = ^2 [pk(qk - qk-i) + H(pk,q k -i)St] 



k=l 



+ Poqo ~PQm - n (p - 1) 



(11) 



and St = t/(M + 1). The last term in this expression 
is specific to the Poisson initial conditions. If the ini- 
tial number of particles is fixed to be no, and therefore 
G(p, 0) = p n ° - the last term is changed to n lnp . The 
same path integral may be derived, of course, using the 
Doi's operator algebra and coherent states. We summa- 
rize this derivation in Appendix ^ The convergency of 
the path integral may be achieved by a proper rotation 
in the complex pk and qk planes. 



FIG. 1: The phase portrait of the binary annihilation process. 
Thick lines represent solution of H(p,q) — 0, fat dot - fixed 
point. Thinner lines represent dynamical trajectories with 
non-zero energy. Line p = 1 gives the mean-field dynamics. 



In what follows we are interested in the semiclassical 
treatment of this path integral. Varying the action with 
respect to pk and qk for k = 0, 1 . . . M, one obtains the 
classical equations of motion (in continuous notations) : 



0H X 2 

q = — 7T- = -*pq ; 
op 

p = ^ = Hp 2 i)? 

and the boundary conditions: 

q(Q) = n ; 
P(t) = P, 



(12a) 
(12b) 



(13a) 
(13b) 



where p and t are the arguments of the generating func- 
tion G(p,t). Notice that while the coordinate is fixed at 
an initial time (past), the momentum is imposed at a fi- 
nite time (future) . These equations admit the integral of 
motion, which we call "energy": E = 0, where 



E = H(p(t),q(t)) = ~{p 2 (t) 



l)q 2 (t). 



(14) 



As a result, the action on a classical trajectory may be 
written as (in continuous notations): 



S[p,q] = Et- qpdt-n o (p(0)-l). 



(15) 



To find the low moments one needs to know G(p, t) in 
the immediate vicinity of p = 1 . In this case the Hamilton 
equations (|12f> with the boundary conditions (|13fl may be 
solved with the mean-field anzatz: 



Pit) 
dq 
~di ' 



1 



dH 
dp 



(16a) 
(16b) 
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The last equation constitutes the mean-field approxima- 
tion for the reaction coordinate, q = n rs (n). The clas- 
sical action, Eq. (|15fl , is obviously nullified on the mean- 
field solution: S[p,q] = 0. This enforces the normaliza- 
tion, Eq. Q (it is straightforward to show that the fluc- 
tuation determinant around the mean-field trajectory is 
unity). In fact, any legitimate Hamiltonian must satisfy 
the condition H(X,q) = to insure normalization. As a 
result, the mean-field solution, p = 1, is bound to have 
zero energy, E = 0. 

However, the assumption that p = 1 is not always a 
legitimate one. The probability of any event other than 
the mean-field prediction is automatically described by 
p(t) = p different from unity. The rare events definitely 
belong to this category. For such cases the mean-field an- 
zatz, Eq. (|I6|I . is not applicable and one must go back to 
the full dynamical system, Eqs. (|12[1 (provided semiclassi- 
cal approximation is justified). For example, let us imag- 
ine doing the contour integral, Eq. , by the stationary 
point method. Approximating, G(p,t) = exp{— S(p, t)}, 
with the classical action, S, one finds for the saddle point 
condition: n = —pdS/dp = p(t)q(t), where we have used 
that on a classical trajectory dS/dp = —q(t). Therefore, 
if one is interested in n which is different from the mean- 
field prediction q(t), one must consider p(t) = p to be 
different from unity. 

In case of the binary annihilation the mean-field pre- 
diction is q(t) = n(t) — no/(l + n^Xt) « (At) -1 for 
1 < (At) -1 <C uq. We are looking for a probability to 
find n 7^ n(t) — (At) -1 particles at time t 3> (Ann) -1 - 
The phase portrait of the dynamical system, Eqs. (fT^|l . 
is plotted on Fig. ^ Dynamical trajectories for a given 
energy, E, are given by q = \J2E\~ 1 / (p 2 — 1). Since 
q(0) = n » 1, one finds p{0) = 1 + 2£/(Ang) » I. 
Substituting this trajectory into Eq. 112b|) and inte- 
grating it between p(0) « 1 and p(t) = p, one finds 
E = — arccos 2 p/ (2 At 2 ). The corresponding classical ac- 
tion, Eq. I|15[l . is given by 
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Fokker-Planck 
Exact solution 
Semiclassical Solution 
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Ratio of the particles number n and average <n> 

FIG. 2: Logarithm of the PDF P n (t) as function of n/n(t) at a 
fixed t. The semiclassical result, Eq. (1181 . - dashed line; the 
exact solution produced by the numerical simulation of the 
Master equation Q - dashed-dotted line; numerical solution 
of the FP equation: P = A[((n 2 - n)P)" + (n 2 P)'] - full line. 



cases the exponent takes the form: 



■ In P„(t) 



n In ZLL ■ n < n , 



— n 

j(n — n) 2 /n; \n — n|<Cn, 
±n /n — n 



(19) 



In 2 ; n>n. 



S(p, t) — i n(t) arccos 2 p . 



(17) 



The logarithm of the PDF is plotted on Fig. versus 
n/n for a fixed n = n(t). The corresponding exponent re- 
sulting from the solution of the Focker-Planck equation 
is shown on the same plot for comparison. The two expo- 
nents coincide for small deviations from the mean-field 
result, \n/n— 1| <C 1. For larger deviations (rare events), 
n/n ~ 0(1), the Focker-Planck results are significantly 
off the correct ones. Finally, the normalization factor 
Af = \Ji/ (Ann) is simply determined by the immediate 
vicinity of the maximum of the distribution, \n — h\ ^ n. 



This action solves the Hamilton-Jacoby equation © and 
is nullified at the mean-field trajectory, p = 1. As 
a result, the generating function is given by G(p, t) k, 
exp{— S(p, t)} with the classical action, Eq. lfP7|l . 

We are now on the position to find the rare events 
statistics: namely we are looking for the probability 
to find n particles after time t, that is P n (t), where n 
is significantly different from the mean field prediction 
n = (At) -1 . To this end one may perform integration, 
required by Eq. J3J, in the stationary point approxima- 
tion to obtain for the probability distribution 



III. BRANCHING AND ANNIHILATION 



Let us consider now a more interesting example of bi- 
nary annihilation with branching. The model consists of 

the two reactions: annihilation A + A and branch- 
ing A 2A. The Master equation is written as: 

d Pn(t) - *[{n + 2)(n + l)P n+2 (t)-n(n-l)P n (t)] 



dt~' iy ~' 2 

+ cr[(n-l)P„_i(t)-nP n (i)] , 



(20) 



P„(t) =A/"exp 



/; | i arccos 2 p s + ^ \np s ) } . ( IN) 



one may check that the corresponding Hamiltonian takes 
the form: 



where p s — p s (n/n) is the solution of the saddle point 
equation: p s (p 2 s — l) -1 ^ 2 arccosp s = n/n. In the limiting 



H(p,q) = ^(p 2 ~l)q 2 -a(p-l)pq. (21) 
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FIG. 3: The phase portrait of the branching annihilation pro- 
cess. Thick lines are lines of zero energy, H (p, q) = 0. Fat 
dots are fixed points. 

As expected, it satisfies the normalization condition, 
H(l,q) = 0. The classical equations of motion are 

q = -Xpq 2 + a(2p- l)q; (22a) 
p = \(p 2 -l)q-a(p-l)p, (22b) 

with the same boundary conditions as in the previous 
example, Eqs. The classically conserved energy is 

E = H(p(t), q(t)). The mean-field anzatz, p(t) = 1, leads 
to the mean-field equation for the reaction coordinate, 
q w (n): 

f t =-X<f + aq. (23) 

This equation possesses two stationary states: the active 
one q = a j\ = n s and the passive one q = 0. Below we 
show that the active state is not actually thermodynam- 
ically stable (in Od system) and in a finite time decays 
into the passive one. 

To proceed with the discussion of the rare events statis- 
tics, we need a phase portrait of the system. It contains 
three lines of zero energy: the mean-field one p = 1; 
the empty system one q = and the non-trivial line 
q = 2n s p/(l +p). These lines determine the topology of 
the phase diagram, Fig.|3| where the arrows show the pos- 
itive time direction. According to the mean-field equa- 
tion 123fl . from any initial state with no ^ 0, the system 
reaches the active state with n s particles during the time 
t w a~ l . Hereafter we assume that n s — a/X ^> 1. We 
shall look for a probability to find n ^ n s particles after 
a time t^a^ 1 . 

Of particular interest, of course, is the probability of 
going to the passive state, namely n = 0, during a large 
time t. According to the definition of the generating func- 
tion, Eq. @, this probability is given by G(0,t). We are 
interested, therefore, in the trajectory which starts at 
some initial coordinate go = (and arbitrary momen- 
tum) and ends at pm = (and arbitrary coordinate) after 
time t. In a long time limit, t — > oo, such trajectory ap- 
proach the lines of zero energy. The system first evolves 



along the mean-field trajectory p = 1 towards the active 
state, q = n s , and then goes along the non-trivial line 
q = 2n s p/ (1 +p) towards the passive state p = q = 0, cf. 
Fig. |3 The action is zero on the mean-field part of the 
evolution, while it is 

o 

/2n s v 
— — dp = n s 2(l-ln2) (24) 
l+p 

l 

along the non-trivial line. 

According to the standard semiclassical description of 
tunnelling 2 -^, to find an escape probability, one has to 
sum up contributions of all classical trajectories with an 
arbitrary number of bounces from (l,n s ) to (0,0) and 
back. Each bounce brings the factor ate~ Sa , where the 
pre-factor reflects the fact the center of the bounce may 
take place at any time without changing the action (zero 
mode). Since the distant (in time) bounces interact with 
each other only exponentially weak, the escape attempts 
are practically uncorrelated. As a result, the probability 
to find an empty system, Po(t) = G(0,t), is 

P (t) = 1 - e-*/ T , (25) 

where the decay time r is given by 

t = <7- 1 cxp{+5 }. (26) 

The semiclassical calculations is valid as long as So 3> 1 
and thus the decay time is much longer than the micro- 
scopic time, t 3> u^ 1 . 

IV. DIFFUSION 

We turn now to the discussion of finite dimensional 
systems. To characterize a microscopic state one need 
to specify number of particles at every site of the lat- 
tice: {ni, . . . , n^}, where N ~ L d is the total number of 
sites. The probability of a given microscopic state may be 
written as P ni ,...,n N (t) and the corresponding generating 
function is 

G( Pl ,..., PN ,t)= ]T P r...p; s P ni „,((). (27) 

ni ,. . .,njv 

Assuming that the reaction rules are purely local (on- 
site), while the motion on the lattice is diffusive, one 
finds that the Hamiltonian takes the form 

H(p 1 ,...,p Nl q 1 ,...,q N ) = ^2 [Ho(pi,qi) + DVpi-Vqi , 

(28) 

where Ho (P> <z) is a zero-dimensional on-site Hamiltonian 
given e.g. by Eqs. Q or J2U; D is a diffusion constant 
and V is the lattice gradient. To shorten notations we 
pass to continuous d-dimensional variable x and intro- 
duce fields p(x) and q(x). The generating function be- 
comes generating functional, G(p(x),t). The latter may 
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be written as a functional integral over canonically conju- 
gated fields p(x, t) and q(x, t), living in d+ 1 dimensional 
space, with the action 



z 

S[p, q] = JdtJ d d x [H (p, q) + DVp -Vq-qp] 



(29) 



The initial term, e.g. the Poisson one: J d d xrio(x)(l — 
p{x, 0)), should also be added to the action. The corre- 
sponding classical equations of motions are: 



q = DV 2 q 



SH 
Sp 



P 



-DX7 2 p + — - 
5q 



(30a) 
(30b) 



These equations are to be solved with the following 
boundary conditions: 



q{x,0) = n (x) : 
p(x,t) = p(x) , 



(31a) 
(31b) 



where uq(x) is an initial space-dependent concentration 
and p(x) is the source field in the generating functional 
G(p(x),t). The mean-field approximation is obtained 
by putting p(x, t) — 1 and is described by the reaction- 
diffusion equation: 



dtq = DV 2 q - 



SH (p, q) 



Sp 



(32) 



that is subject of numerous studies. 

Equations l|30|l admit the integral of motion: E = 
J d d x [H (p, q) + D S7pVq] . In some cases (see below) an 
additional infinite sequence of integrals of motion may 
be found, making the classical problem, Eqs. (J20J), an- 
alytically solvable. In a general case, these equations 
must be solved numerically. We notice, however, that 
such numerical problem is orders of magnitude simpler 
than numerical solution of the Master and even FP equa- 
tions, or direct modelling of the stochastic system. Below 
we discuss a fast, efficient algorithm for numerical solu- 
tion of Eqs. (|30|1 with the boundary conditions Eqs. <|31|) . 
Moreover, a lot of insight may be gained by investigating 
the phase portrait of the zero-dimensional Hamiltonian, 
Ho(p,q), which allows to make some semi-quantitative 
predictions without numerical solution. 

To illustrate how the method works we consider the 
branching annihilation problem of section lTTTl (Hq is given 
by Eq. 1121(1 ) on a compact d-dimensional cluster - the 
"refuge"—, denoted as TZ. Outside of the refuge, there is 
a very high mortality, A — > 0, rate which is eventually 
taken to infinity. This dictates the boundary condition 



q{dK,t) =0, 



(33) 



where &TZ is the boundary of the cluster TZ. It is conve- 
nient to pass to the dimensionless time at — > t and coor- 
dinates x/£ — > x, where £ = \jD/a. We also introduce 



the rescaled fields q(x,t) = n s ip(x,t) (where n s = er/A) 
and p{x,t) = 1 — 0(x, t). In these notations the semi- 
classical equations, Eq. I(3()|l . take the symmetric form: 



d t ip — V 2 ip + tp — ip 2 + (pip 2 — 2iptp , (34a) 
-d t <p = \7 2 tp + tp-tp 2 + ip>p 2 - 2pip , (34b) 

Consider first the mean-field (ip — 0) evolution, de- 
scribed by the equation 



d t p = VV + p - p 2 



(35) 



subject to the boundary condition ip(dlZ,t) — 0. For 
the small concentrations, p -C 1, the last term may be 
omitted and the solution takes the form: 



<p(x,t) = Ya n e( 1 - x ^ t Y n (x) 



n=0 



(36) 



where Y n (x) are normalized eigenfunctions of the Laplace 
operator in the region TZ with zero boundary conditions 
and eigenvalues — A„ < 0; coefficients a n depend on an 
initial condition. Therefore, if the smallest eigenvalue, 
Ao, is larger than unity (the cluster is small enough), 
any initial distribution evolves towards the empty sys- 
tem. The characteristic lifetime of the system is thus 



^(Ao-l)- 



A > 1. 



(37) 



If Ao < 1 (the cluster is larger than some critical size), 
the mean-field evolution, Eq. I|35l) , leads to a stable non- 
vanishing concentration po (r) , which is given by the solu- 
tion of the equation ~S/ 2 ipo + po — p\ = with zero bound- 
ary conditions. It is clear, however, that such solution is 
actually a meta-stable state of the system. Namely, af- 
ter a long enough time the system will find itself in the 
empty (passive) state. Our task is to find the system's 
lifetime, r, for the meta-stable case, Ao < 1. According 
to our previous discussions the lifetime is expected to be 
exponentially long 



-i e s d 



Ao < 1, 



(38) 



where Sd is the action along the semiclassical trajectory, 
that solves Eqs. (|34a|l . I|34b|) with the initial condition 
ip(x,0) = <po(x) and the final condition p(x,t e ) = 1. 
The extinction time, t e , is to be sent to infinity. Indeed 
dSd/dt e — E(t e ) < and thus the longer the extinction 
time - the smaller the action. In practice, however, the 
action almost saturates at modest values of t e . 

In general the problem cannot be solved analytically 
and one needs to resort to numerical approaches. The 
following iteration scheme rapidly converges to the de- 
sired solution: one first fixes momenta to be p>i(x, t) = 1 
at any time and solves Eq. (|34a|) with the initial condi- 
tion p{x 1 0) = (po(x) by forward iteration from t = to 
t = t e . The result of this procedure, pi(x, i), is kept fixed 
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Radius, R in \ units 

FIG. 4: The semiclassical action (in units of £ d n s ) for the 
extinction of a one dimensional cluster is shown as function 
of cluster's radius 7? (in units of f) - full line. Large radius 
approximation, Eq. 1391 . is shown by the dashed line; near- 
critical, 7r/2 < R, approximation, Eq. I|450 - dashed-dotted 
line. 

during the next step, that is solution of Eq. I|34b(l with 
the condition if(x, t e ) = 1 by the backward iteration from 
t = t e to t = 0. This way one finds ^(x, t), which is kept 
fixed while the next approximation Lp2(x,t) is obtained 
by the forward iteration of Eq. I|34a[l . Repeating suc- 
cessively forward and backward iterations the algorithm 
rapidly converges to the required solution. The action 
Sd = Sd{t e ) is then calculated according to Eq. (|2*9")) . 
Finally one has to check that Sd(t e ) does not decrease 
significantly upon increasing t e . 

The action, Sd, for one-dimensional cluster of size 2R 
is plotted on Fig. 21 as function of R. The critical ra- 
dius R c = tt/2 (in units of £) is found from the condi- 
tion Aq = 1. For R < R c , there is no meta-stable state 
and thus Sd — 0, while the cluster lifetime is given by 
Eq. For R > R c the lifetime is given by Eq. (|38[l . 

with the numerically calculated Sd plotted on Fig. 0] 
The asymptotic behavior of the action Sd for R ^> R c 
(A < 1) and R - R c < R c (1 - A < 1) may be readily 
found analytically. 

For R ;§> R c (Xq <C 1) the concentration throughout 
the bulk of the cluster is practically uniform apart from 
a surface layer of thickness £. One may therefore apply 
the results of the zero dimensional problem, Eq. H24|) . to 
find 

S d = 2(1 - \n2)n4 d (V - cS) , (39) 

where V and S are cluster's dimensionless volume and 
surface area correspondingly and c is a numerical con- 
stant, which we shall not evaluate here. For the Id case 
the corresponding line is plotted on Fig.^Jby the dashed 
line. 

We turn finally to the clusters that are only slightly 
larger than the critical size: e e 1 - Aq C 1 (e.g for 



a spherical cluster with the radius R one finds e = 1 — 
(R c /R) 2 , where R c is a critical radius, found from Ao = 
1). In this case only the zeroth eigenfunction Yo(x) is the 
unstable direction of the linearized mean-field equation. 
The full non-linear mean-field equation l|35|l possesses, 
therefore, the stable solution ipo(x), that is expected to 
be of order e. One may thus look for this solution in the 
following form: 

Mr) = e(r)Y (x) + e<p 1 (x)) , (40) 

where (fi(x) is orthogonal to Yq(x). One can now sub- 
stitute this trial solution in Eq. I)35|l . keeping only the 
leading (second) order in e, and project on Yq, using its 
orthogonality to (pi . As a result the coefficient 77 is found 
to be: 

rf x = Jd d rY^{r). (41) 
n 

The meta-stable solution of equations $5<^i in the leading 
order in e is therefore: ifio(x) = er]Yo(x) and (f>o(x) = 0. 
To find the optimal escape trajectory, let us parameterize 
deviations from this meta-stable state as 

00 

<p(x,t) = er)Y (x) + J2an(t)Y n (x); (42a) 

00 

ip(x,t) = J2Pn(t)Y n (x), (42b) 

n=0 

where a n (t) and f} n {t) are assumed to be small. One 
can now substitute these deviations into the dynamical 
equations (|34|l and linearize them with respect to a„, (3 n . 
It is straightforward to see then that in the leading order 
in e only ao and /3n should be retained. They evolve 
according to: 

!(£)=■( 1 2 )(?:)+ <*■>■ <«) 

The matrix on the right hand side has two eigenvec- 
tors, (1,0) and (1,-1) with the eigenvalues —1 and 1 
correspondingly. The first one describes deviation in the 
mean-field direction, ip = 0, and leads to the restoring 
force back to the meta-stable state. The second eigen- 
vector gives the most unstable direction, that describes 
the way the system escapes towards the empty state. The 
corresponding trajectory on the (^?, ip) plane is plotted on 
Fig- El for the center point of the Id cluster, x — 0. Dif- 
ferent lines correspond to a few values of t e . For t e — » 00 
the energy, E, approaches zero and the trajectory ap- 
proaches the (1,-1) direction that leads from the meta- 
stable point (0,e?7Yo) to a symmetric meta-stable point 
(erylojO). (The existence of the latter follows directly 
from the symmetry of equations (|34|l .') For e 1 the 
small deviation analysis describe the entire transition be- 
tween the two meta-stable points, that takes place, there- 
fore, along the straight line: 

<p{x,t) = er)Y (x) - tp(x,t) , (44) 
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FIG. 6: The phase portrait of the run-away process, Eq. 1461 , 
Thick lines represent solution of H(p,q) — 0, fat dots - fixed 
points. 8 2 = 1/2. 



FIG. 5: Trajectories on the (<p, <p) plane is shown for the cen- 
ter point of the Id cluster, x = 0. The horizontal axes is ip 
variable and the vertical axes is cp. Different curves are distin- 
guished by their escape time, marked in dimensionless units 
at e . er)Y o (0) = 0.3. 



on the ((/?, if) plane. The farther evolution takes place 
along ip — direction. As a result, in the limit t e — > oo, 
and therefore E — > 0, the semiclassical escape action is 
given by the area of the straight triangle with the hight 
er)Yo(x) integrated over the cluster, cf. Eq. J2HJ), 



S d 



TZ TZ 



1 



x hpdip = - {r)e) 2 n s t; d 



(45) 

For the Id cluster (e = 1 - (tt/2R) 2 and rj 2 = 9?r 3 /128) 
Eq. I|45l) is shown on Fig. EJby a dashed-doted line. (For 
a circular cluster in 2d R c = 2.4 and rj 2 = 9.4, while for 
a spherical 3d cluster R c — tt and rj 2 — 51.7.) One may 
observe that the large and small cluster asymptotic re- 
sults, Eqs. (|59")l and illST) correspondingly, provide a rea- 
sonable approximation for the exact numerical calcula- 
tion of the semiclassical action, Sd- Finally, the prob- 
ability of the system staying in the meta-stable state 
is P(t) = exp{— t/r}, where the lifetime r is given by 
Eq. J3BJ). 



tion for the zero-dimensional system has the form: 



dP n 

dt 



A 



(n-l)(n-2) n(n-l) D 

— r" n —l "*> 



2 2 
a \{n + l)P n+1 - nP n ] + n [P„_ x - P n ] .(46) 



The corresponding zero-dimensional Hamiltonian is: 
H (p, q) = ^(p 2 - f)q 2 + a(p - l)q + M (l - p) (47) 

and the classical equations of motions are: 

3 



q = -X(p- 7^P 2 )q 2 -aq + fJ,] 
p = X(p 2 -p 3 )q + a(p- 1). 



(48a) 
(48b) 



As always, the mean-field equation of motion for the re- 
action coordinate q w (n) is obtained by the anzatz p = 1 
and takes the form 



dq A _ 2 



(49) 



According to the mean-field equation there are two qual- 
itatively different scenarios of the system's evolution. 
They are distinguished by the parameter 



2Xfj, 

,t2 



(50) 



V. RUN-AWAY SYSTEMS 

In this section we consider a qualitatively different sys- 
tem that exhibits a run-away behavior, characterized by 
unlimited proliferation of the number of particles. The 
simplest example is given by the population dynamics 
model consisting of three reactions: binary reproduction, 
death and emigration, characterized by probabilities A, a 
and \i correspondingly. The schematic way to write it is: 



A + A 



3A; A 



and 



A. The Master equa- 



If 5 2 < 0, the r.h.s. of Eq. (|4"9"|l is strictly positive and the 
reaction coordinate always grows to infinity. This is the 
scenario, where the population proliferates indefinitely. 
Alternatively, for 8 2 > the system possess two station- 
ary concentrations: n T = n s (l =p S), where n s = cr/X. 
The point q — n_ is the stable one, while q = n+ is un- 
stable. In this case (the only one we discuss hereafter), 
the mean-field predicts that for the range of initial con- 
centrations < tiq < n + the system evolves towards the 
stable population n_ . If the initial concentration exceeds 
n+ - the system runs away and the population diverges. 
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If one goes beyond the mean-field treatment, how- 
ever, one realizes that the state n_ is actually a meta- 
stable one. To see this fact and calculate the life-time 
of the meta-stable state, it is convenient to draw the 
phase portrait, Fig. |(|] It has two lines of zero energy: 
the mean-field one, p = 1, and the non-trivial line 
Xp 2 q 2 /2 — aq + p = 0. These two lines intersect at the 
mean-field stable points p = 1 ; q = and determine 
the topology of the phase diagram. It is clear from the 
phase portrait that the point p — 1 ; q — n_ is not 
stable once motion with p ^ 1 (non-mean-field) is al- 
lowed. More precisely, there is a non-mean-field path 
that brings the system from the point q — n_ to the 
point q = n+ . Once the point q = n+ is reached, the sys- 
tem may continue to evolve according to the mean-field 
towards indefinite population grow. Repeating calcula- 
tions, similar in spirit to calculations of the decay-time in 
section ITTT1 one finds for the life-time of the meta-stable 
state, q = n_: 



exp{+S } , 



(51) 



where Sq is the classical action along the non-trivial line 
of zero energy between points (1, rt_) and (1, n+). Calcu- 
lating the integral, one finds So = f( n +) ~ f( n -)i where 
f(x) = x — y/%fi/X arctan(xyA/2/I). 

Two limiting cases are of particular interest: (i) the 
"near critical" system, < S 2 <C 1 and (ii) the sys- 
tem with almost no immigration, fi — > + ; 5 — > 1 _ . 
In the former case the two mean-field stationary points 
approach each other, making the escape from the meta- 
stable state relatively easy. Expanding the /-function up 
to the third (!) order, one finds Sq — 2n s S 3 /3 -C n s . 
As expected, the action is small and correspondingly 
the life-time is short (notice that the quasi-classical pic- 
ture applies as long as So > 1). In the latter case 
the two mean-field stationary points tend to n_ — > 0, 
and n+ — > 2n s . If the immigration is absent, p = 0, 
the mean-field stable point, n_ = 0, coincides with the 
empty state of the system. The empty state is abso- 
lutely stable since no fluctuations are possible. Naively 
one may expect that in this limit the life-time of the 
meta-stable state (and thus So) diverges. This is not the 
case, however. The calculation shows: Sq — > 2n s . As a 
result, even negligibly small probability of immigration, 
H, leads to a finite probability of unlimited population 
expansion. (Strictly speaking, one also needs to show 
that the pre-exponential factor does not go to zero once 
M-0.) 

We consider now a finite-dimensional generalization of 
this population dynamics model. The physics of the phe- 
nomena, discussed here, is as follows: if a critically large 
cluster "tunnels" into the run-away state, both diffusion 
and reaction dynamics work to expand the cluster and 
flip the entire system into the run-away mode. The situ- 
ation is similar to nucleation of a critical domain in the 
super-cooled state of a system close to a first order phase 
transition. To simplify the algebra we shall consider only 
the case of the " near critical" system, < 8 2 <C 1 , where 



the apparatus turns out to be rather similar to that of 
the theory of the first order phase transitions. 

As discussed above, the finite-dimensional generaliza- 
tion of the Hamiltonian is H[p,q] = fd d x[Ho(p,q) + 
DVpVq]. For 5 <C 1, it is convenient to make a change of 
variables (p, q) — > {p,p), as p = l + ip and q = n s (l + ip), 
where ip ~ S, while p ~ S 2 . Substituting it into the reac- 
tion part of the Hamiltonian, Eq. 147|) . and keeping terms 
up to S 4 , one obtains Ho(p, ip) = an s [tp(5 2 — p 2 )/2 — tp 2 ]. 
As a result, the d-dimensional action, Eq.JSHJ, for the 
conjugated fields tp(x,t) and p(x,t) takes the form: 



S = n s ^ d dt 



d d 3 



tp ( ip — V (p + 



(52) 

where we have introduced the dimensio nless time at — > t 
and coordinate x/£ — > x, where £ = y/D/cr. The func- 
tional integration over the field tp should be understood as 
running along the imaginary axis. The field theory with 
the action Eq. Il52[l may be considered as Martin-Sigia- 
Rose^l representation of the following Langevin equation: 



dip 



V 2 



tp ■ 



dV 

dp 



(53) 



where t) is a Gaussian noise with the correlator 

(C(x, t)C(x', t')) = ^ 5(x - x')5{t - t') (54) 

and the potential is V(p) = —p 3 /6 + S 2 p/2 . This 
potential has a meta-stable minimum at ip = —S and 
an unstable maximum at tp = 8. The barrier hight is 
V{6) - V(-6) = 26 3 /3 and therefore the life-time of the 
zero dimensional (d — 0) system is expected to be given 
by the activation exponent (with (n s ^ d )^ 1 playing the 
role of temperature) ~ exp{rt s 2<5 3 /3}, in agreement with 
Eq (£>D. 

To discuss the life-time of the finite-dimensional sys- 
tem we shall not use the Langevin approach, but rather 
return to the action, Eq. (|52|l . and write down the clas- 
sical equations of motion: 



dtp = VV 



dV 

dp 



2p; 



d t p> 



-V 2 ip + ip 



d 2 V 
dp 2 



(55a) 
(55b) 



The energy density, corresponding to these two equa- 
tions, is defined as E(x,t) — —ip(W 2 p — dV/dp + tp). 
The global energy, E = Jd d x E(x,t), is, of course, con- 
served. However, in the present case if E(x, 0) = 
it keeps holding locally at any time: E(x,t) = 0. In- 
deed, the energy density vanishes if either tp = 0, or 
tp = — V 2 p + dV/ dp = 2p — dtp and thus tp = dtp, where 
we have employed Eq. I|55a|) . It is easy to check that in 
both cases Eq. (|55b(l is satisfied automatically. There- 
fore the evolution with zero energy density is described 
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by either dtf — V 2 ip — dV/d(p, that is the mean-field 
equation, or by d t cp = —\7 2 ip + dV/d<p, that gives the 
motion along the non-trivial line of zero energy. 

Notice that the last equation happens to be the time- 
reversed version of the mean-field^. If one starts, thus, 
from the stationary solution ip = — 8 and perturbs it in- 
finitesimally - the perturbation grows until it reaches the 
stable configuration, satisfying 

VV-^=0. (56) 

The critical domain is given therefore by a localized so- 
lution of Eq. I|56[). Since the energy along the nucle- 
ation dynamics is zero, the action to nucleate the crit- 
ical domain is given by S d — n s t; d fd d x fdtipdttp = 
n s t; d J d d x Jdt (—V 2 p+dV/ dip) dt<p. Performing the time 
integration in this expression, one finds for the action 

S d = n4 d J d d x Q(V^) 2 + V(tp d ) - V(-5)) , (57) 

where ip d = ipd(x) is a stationary localized solution of 
Eq. (|5(j|) , that is an extremum of the functional (|57|l . As 
a result, the problem of the dynamical escape from the 
meta-stable configuration is reduced to the static Landau 
theory of the first order transitions. As far as we know, 
such reduction is not a general statement, but rather is 
a consequence of the assumption 8 <C 1 and the result- 
ing local energy conservation, E(x, t) =0. In a general 
situation, one still has to solve a considerably more com- 
plicated problem of dynamic equations l|48|) for <p(x,t) 
and (p(x, t). 

^From the scaling analysis of Eq. i|56|) one finds that 
ip ~ 8 in the core of the critical domain. Employing this 
fact, one finds that the characteristic spatial scale of the 
domain is given by 8~ x l 2 3> 1 (distance is measured in 
units of £ = yj D/a ) . Therefore the action cost to create 
the critical domain is 

S d = c d n s (?p- d > 2 , (58) 

where Cd is a numerical factor of order of one: Co = 
2/3 ; ci = 24/5. This result suggests that for d > 6 the 
state with finite population density n — n s (l — 8) is sta- 
ble, while for d < 6 the state is meta-stable. The concen- 
tration of critical domains is given by £~ d exp{— Sd} and 
the typical distance between them is £ exp{S d /d}. They 
grow diffusively until the entire system is flipped over to 
the run-away state in time r ~ er -1 exp{2S d /d}. The 
semiclassical calculation is applicable as long as S d > 1 
and therefore 8 is not too small. For very small 8 the 
escape is driven by the fluctuations rather than the semi- 
classical dynamics. 

VI. CONCLUSIONS 

The examples, considered above, are meant to illus- 
trate the general technique to calculate probability of 



rare events in reaction-diffusion systems. The technique 
is based on the existence of the many-body "quantum" 
Hamiltonian, which fully encodes the microscopic Master 
equation. The very same Hamiltonian in its second quan- 
tized representation serves as a starting point for field- 
theoretical treatments of dynamic phase transitions in 
the reaction-diffusion system 11 ! 12 ! 1 ?! 14 ! 1 ^ 11 ?! 1 ?! 1 ?! 1 ?. For 
our present purposes we have deliberately chosen to work 
with systems that are away from a possible continuous 
phase transition point. Namely, we focus on the parts 
of the phase diagram, where the mean-field considera- 
tions suggest a non-vanishing population of particles (or 
at least transiently non- vanishing population). In such 
cases the "quantum" fluctuations are small and one may 
treat the underlying "quantum" dynamics in the semi- 
classical way. 

We stress that the semiclassical treatment is not equiv- 
alent to the mean-field one. The latter requires a very 
special assumption about dynamics of the canonical mo- 
menta: p(x,t) = 1. This assumption may be justified 
as long as one is interested in a typical system's behav- 
ior (even this is not guaranteed if the system possess 
meta-stable states, as in our last example). In such cases 
the problem is reduced to a partial differential equation 
for the reaction coordinates, q(x,t), only. However, if 
questions about atypical, rare events are asked, - the 
mean-field assumption, p(x,t) = 1, must be abandoned. 
As a result, one has to deal with the canonical pair of 
the Hamilton equations for reaction coordinates, q(x,t), 
and momenta, p(x,t). The degree of deviation from the 
mean-field line, p — 1, is specified (through proper ini- 
tial and finite boundary conditions) by the concrete sort 
of the rare event of interest. Finally, the probability of 
the rare event is proportional to the exponentiated action 
along the classical trajectory, satisfying specified bound- 
ary conditions. 

We found it especially useful to work with the phase 
portrait of the corresponding dynamical system on the 
(p, q) plane. The emerging structures are pretty intuitive 
and can tell a great deal about qualitative behavior of 
the system even before any calculations. The Hamilto- 
nians underlying the Master equations of reaction sys- 
tems are typically not of the type traditionally consid- 
ered in the theory of dynamical systems. For exam- 
ple, they usually can not be cast into the familiar form 
H(p, q) = p 2 /2 + A(q)p + V(q). On the other hand, they 
posses some universal features, such as H(l,q) = 0, or, 
if there is an empty absorbing state, H(p, 0) = 0, etc. 
These features dictate a specific topology of the phase 
portrait. It would be extremely interesting to explore 
this class of Hamiltonians from the point of view of math- 
ematical theory of dynamical systems^. A question of 
particular interest is a possible exact integrability of re- 
sulting Hamiltonian equations (especially in d = 1)24. 

There are number of issues, that are not addressed in 
the present paper and require further investigation. Let 
us mention some of them, (i) Throughout the paper we 
have discussed the rare events probability with the expo- 
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nential accuracy. In some cases this is not enough and 
one wants to know the pre-exponential factor rather pre- 
cisely. This requires calculation of the fluctuation deter- 
minant on top of the non-trivial classical trajectory. This 
task is relatively straightforward for the d = systems, 
where it may be addressed by writing down "quantum" 
corrections to the Hamilton- Jacoby equation and treat- 
ing them iteratively (in the way it is usually done in the 
single-particle WKB method) . For extended systems the 
task is reduced to the spectral problem of a certain ma- 
trix differential operator. At present we are not aware of 
a general recipe to solve it. One may show, however, that 
on any mean-field trajectory, p{x,t) — 1, the fluctuation 
determinant is equal to unity. The simplest way of doing 
it is to use the discrete representation of the functional 
integral, Eq. (|10|) . and notice that the quadratic fluctua- 
tion matrix has a triangular structure with unities on the 
main diagonal (and, hence, unit determinant). Unfortu- 
nately, this is not the case away from the mean-field, 

(ii) We have restricted ourselves to the systems with 
a single sort of spices only. It is straightforward to 
generalize the technique to any number of spices, K. 
The difficulty is that the phase portrait becomes a 2K- 
dimensional construction, which is not easy to visual- 
ize. Correspondingly, the mean-field line becomes K— 
dimensional hyperplane. Moreover some qualitatively 
new physics may arise such as stable oscillatory limiting 
cycles on the mean-field hyperplane. A paradigm of such 

behavior is a Lotka-Volterra^ system: A + B — > 2A; 
A and B — — > 2B. An example of rare event may 
be an "escape" from the periodic limiting cycle on the 
A — B mean-field plane into the empty state in a finite 
size system. Finding an optimal "reaction path" for such 
escape is not a straightforward matter, however. 

(iii) We have not treated long-range interactions 
and (local or non-local) constraints. The simplest 
(" fermionic" ) constraint is that of a maximum single oc- 
cupancy of each lattice site. It was shown recently that 
such constraint may be incorporated into the "bosonic" 
formulation 36 , leading to a new class of the interesting 
Hamiltonians. Studding rare event statistics for such 
hard-core particles (by studding classical dynamics of the 
corresponding Hamiltonians) is a very interesting sub- 
ject. 

(iv) There is a close resemblance between the for- 
malism presented here for essentially classical systems 
and the Keldysh technique for non-equilibrium quantum 
statistical The semiclassical solutions with 1, con- 
sidered here, correspond to saddle point configurations 
of the Keldysh action with a different behavior on the 
forward and backward branches of the time contour. Al- 
though examples of such saddle points were considered in 
the literatur o 38 ! 39 , it would be interesting to learn more 
about possible applications of the present technique for 
true quantum problems. 

We are grateful to A. Elgart, Y. Gefen, A. Lopatin and 
K. Matveev for useful conversations. A. K. is A. P. Sloan 
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APPENDIX A: OPERATOR TECHNIQUE 

We give here a brief account of the operator 
techniqueS&SiSSiffl for completeness. Define the ket- 
vector \n) as the microscopic state with n particles. Let 
us also define vector 

oo 

|*(t))=£p n (t)|n>. (Al) 

Notice, that the weight, P n , is probability rather than 
the amplitude. It is convenient to introduce the creation 
and annihilation operators with the following properties: 

a f |n) = |n+l); (A2a) 
a\n) = n\n - 1) . (A2b) 

As a byproduct, one has a\0) = 0. One may immediately 
check that such operators are "bosonic" : 

[a, a f ] = 1 . (A3) 

As for any pair of operators satisfying Eq. 1|A3(I one may 
prove the identity 

e Q /(a,« t ) =/(a,a t + l)e°, (A4) 

where / is an arbitrary operator-value function. In these 
notations the whole set of the Master equations may be 
recast into single " imaginary time" Schrodinger equation 

^|*(i)> = -#!*(«)>, (A5) 

where H is the "Hamiltonian" operator. One may check 
that the Hamiltonian of the binary annihilation process, 
Eq. QJ, has the form 

H=^f-l)a\ (A6) 

where the first term in brackets on the r.h.s. is the "out" 
and the second one is the "in" term. 

One may solve formally the Schrodinger equation and 
write |*(t)) = exp{-i?(a t , a)t}|*(0)). An initial state, 
|*(0)), is specified as e.g. |<f(0)} = e -"o(a + -i) | ) for the 
Poisson initial distribution, or |^(0)) = (a^) |0) for the 
fixed particle number. The generating function Eq. (J2J 
is given by 

G(p,t) = (0\e pa e - 6 ^' a ^(0)). (A7) 

The normalization, (7(1, t = 0) = 1, is guaranteed by the 
identity (0|e a |n) = 1 for any n (this fact may be checked 
using Eq. $M$) and the fact J2 n P n {0) = 1- The nor- 
malization is kept intact at any time if {0\e a H (a ! , a) = 0. 
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Since the coherent state (0|e a is an eigenstate of the cre- 
ation operator, (0|e a a t = (0|e a , one arrives at the con- 
clusion that any legitimate Hamiltonian must obey 

H{a) = l,o) = 0. (A8) 

E.g. the Hamiltonian of the binary annihilation, 
Eq. HA6JI . indeed satisfies this condition. 



One may employ now the standard bosonic coher- 
ent state technique to write the generating function, 
Eq. I|A7(I . as the functional integral. The result coin- 
cides identically with the Eq. (|l(Jf> of the main text. One 
notices, thus, the formal correspondence between the op- 
erators a) and a and operators p and q correspondingly. 
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